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1 Introduction 

The flow and handling of granular materials is of major importance to many 
industries. Yet despite efforts over several decades, the modeling of such flows 
has achieved only modest success. Dense gravity-driven flows in hoppers have 
been often modeled as elastic-plastic continua, for example. In this picture, the 
granular material flows as a plastic with a frictional yield condition, and de- 
forms as an elastic solid otherwise. This model has been used to analyze mass 
flows, where the material is flowing throughout the hopper, but has failed to 
provide adequate agreement with experiment in the prediction of quantities such 
as discharge rate, for example Despite such shortcomings, Jenike's || con- 
struction of steady-state incompressible rigid-plastic radial solutions in hoppers 
with simple geometries has been of considerable importance in hopper design 
[ p"o[ . These are solutions for quasistatic flow (inertial effects are neglected) 
where grains travel radially in conical or wedge-shaped hoppers. Only recently 
have numerical solutions of steady-state quasistatic flows in more complicated 
geometries been produced || . 

However, it is now clear that there are serious mathematical difficulties with 
the equations for time-dependent incompressible rigid-plastic flow (IRPF). In 
many instances, the equations for such flows have been shown to be ill-posed 
i.e. they possess instabilities with arbitrarily short wavelengths (Schaeffer Jl2[ , 
Valanis and Peters [|l8)). Hence, it is problematic to interpret the steady-state 
rigid-plastic flow equations as long-time solutions of the time-dependent equa- 
tions. Additionally, both the steady-state and time-dependent equations have 
physical shortcomings. There are several features of hopper flows where the 
particle size is important, indicating that it may not be appropriate to model 
such dense granular flows as continua. 
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Figure 1: Schematic illustration of shear banding. In granular materials the 
band is typically several grain diameters (d) thick. 

One flow feature that involves the particle size is that of shear-banding. 
Shear-banding is an extremely interesting phenomena that occurs in plastically 
deforming materials, including granular materials. At high strains the deforma- 
tion of the material becomes localized in thin bands of shearing. This leads to a 
jump in velocity across the band as illustrated in figure [l| In granular flows, the 
thickness of this band is typically only several grains thick, yet the length can be 
comparable to the size of the container. Thus shear-banding in granular flows 
seems to be a flow feature that spans both the kinetic and the hydrodynamic 
regimes. In this sense, shear-banding in granular materials differs, for example, 
from the discontinuities encountered in compressible fluid flows where shocks 
are smoothed by hydrodynamic effects (e.g. viscosity). The full theory of gran- 
ular flow ought to predict a shear band thickness of zero in the hydrodynamic 
limit. 

It is probable that a complete description of shear-banding in granular flows 
is not within the scope of a continuum theory. Nonetheless, attempts have 
been made to describe shear-banding in continuum theories in some average 
sense (Q). Indeed, continuum plastic flow theories have been found to ex- 
hibit shear-banding type behavior and the development of shear-bands seems 
to be closely related to the development of ill-posedness in such flows jl4j|. 
One approach to deal with ill-posedness is to model the granular material as a 
Cosserat continua. Cosserat-type continuum theories attempt to model materi- 
als with internal structure by the inclusion of extra terms motivated by physics 
at the granular scale. These extra terms can damp the instabilities that lead 
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to ill-posedness and can be used to predict shear-band thickness (Muhlhaus 
[pi). Including such damping terms which act at a granular length-scale and 
'thickening' the shear-bands is a practical approach to studying granular flows 
using continuum equations. Indeed, such an approach is analogous to including 
Reynolds stresses in the Navier-Stokes equations to model turbulence. 

The purpose of this paper is to study a set of equations for granular flow, 
which are a simple generalization of the IRPF equations, and which display 
shear-banding instabilities similar to those that arise in the Critical State theory 
(CST) of granular flow |r|. These equations appear to be more regular than the 
full IRPF equations and are somewhat simpler than CST-motivated equations. 
We do not explicitly introduce a granular length scale here; we wish to study 
the development of shear-banding from a hydrodynamic viewpoint. In the first 
section we state the usual equations for IRPF. In the next section , we review 
the linear stability analysis of these equations and show that they are linearly 
ill-posed in two dimensions. In section four, we introduce a new set of equations 
which are a singularly perturbed form of the IRPF equations, and in the next 
section we present numerical solutions of these equations which demonstrate 
shear-banding. Finally, we look at the linear stability of these equations to 
provide a qualitative description of the shear-banding process. 



2 Equations for an incompressible rigid-plastic 
flow 

We consider the flow under gravity of an incompressible granular material in a 
wedge-shaped hopper under plane strain. We model the flow of this material 
as a continuum rigid-plastic flow. The equations for such a flow consist of the 
incompressibility condition, 

V;^ = G (1) 
where u l are the components of the velocity field, and the momentum equations: 

pi-^+u^ju*) +ViO» =pg\ (2) 

where Oij is the symmetric stress tensor, and g l is the acceleration due to gravity. 
Note that we define o to be positive when forces are compressive. 

Plastic deformation is assumed to occur everywhere in the hopper i.e. the 
material is at plastic yield throughout the hopper. We will use a flow rule 
based on a Drucker-Prager type yield condition . Specifically, in terms of the 
principal stresses Oi, this condition is written as 

{a x - of + (o 2 - of + (o 3 - of < k 2 o 2 (3) 

where k — v^sin^, ip is the internal angle of friction of the material and 
o = \ {o\ + (72 + 03) is the average trace of the stress tensor (we will refer to o 
as the average stress). If this inequality is satisfied exactly then the material is 
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deforming plastically. Under plane strain a 2 = a = h(ai + 03), so in this case 
we may consider a strictly two-dimensional yield condition: 

(0-1 - a) 2 + (fj 3 - a) 2 = 2a 2 sin 2 (p (4) 

We now assume a flow rule of the form 

<7y = (T(5y + /iV^, (5) 

where = (VjUj + VjUj)/2 and /it is some, as yet unspecified, scalar function 
of the normal stress and strain rates. If we compare the flow rule (^) to the yield 
condition (|^), then we see such a flow will satisfy the yield condition exactly if 
we choose the function fi to be 

^ = 777717- ( 6 ) 



where ||V|| = JV~W. 



Equations (Q) , , (||) and (Q) form a closed system for incompressible rigid- 
plastic flow in plane strain. For granular flows in hoppers, these equations are 
only valid for so-called mass flows where the material is flowing throughout the 
hopper. When the hopper is not sufficiently steep funnel flows can develop where 
material flows down a central funnel leaving a stagnant region adjacent to the 
walls. Indeed, radial solutions have been used to study the transition between 
mass and funnel flow which is thought to occur when the rigid-plastic equations 
become singular as the rate of deformation vanishes We will restrict our 
attention to mass flows where rigid-plastic flow occurs everywhere in a given 
domain. 

On each wall of a hopper, the direction of flow is tangent to the wall and 
the stresses must satisfy a Coulomb friction condition. The Coulomb boundary 
condition in two dimensions can be written as follows: 

(Jijrft 3 = tantpu, OijTi l n > , (7) 

where n is a outward-pointing normal vector to the wall, t is a tangent vector 
to the wall pointing in the opposite direction to that of the flow and ip w is the 
angle of wall friction. 

Combining equations ([!]), (j|), @ and (^) together we obtain the equations: 

— +u^ 3 u l = g l -V 3 { P 5v -kpA«), (8) 



where 



V lU * = 0, (9) 



P = cr/P, 



These equations are strictly hyperbolic in the case of plane strain (equation (El)). 
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3 Ill-posedness of the incompressible rigid-plastic 
flow equations 

The equations (||) and have been shown to be linearly ill-posed in certain 
geometries and for certain parameter values Jl^| . Specifically, the linearized 
equations of motion that describe the propagation of a small disturbance in 
the flow, possess unstable plane-wave solutions in the short wavelength limit. 
It has been suggested that these instabilities are related to the formation of 
shear bands in the flowing granular material. In this section, we will consider 
the linearized equations of motion for a plane-wave disturbances to look for 
ill-posedness. This ill-posedness in the IRPF equations was originally noted by 
Schaeffcr |l^| for fully three-dimensional flows. Here we will work only in 2D, 
but with the full linearized equations of motion for a small disturbance. This 
will serve to provide a comparison with work in the following sections. 

We will begin by writing the two-dimensional rigid-plastic equations and 
^|) in non-dimensional form as follows: 

u 1 — u 1 /uq, p — p/gL, x l =x l /L, i — tu /L, (10) 

where uq is some characteristic velocity and L is some characteristic length-scale 
of the problem. The equations for the rigid-plastic flow then become 



du l 
~dt 



^(ffVff-Vi^tf-fcMy]), (11) 
V.m 1 = 0. (12) 



where Fr — uo/^/gL is the Froude number. At this stage will drop the " notation 
and assume that we are dealing with dimensionless quantities unless otherwise 
specified. 

The linearized equations of motion for a small disturbance (du, Sp) propa- 
gating on a smooth background solution (u,p) to equations ( |ll|) and ( |l2l ) can 
be shown to be: 

-^- + u 3 Vj<5u 1 + *u , V J V = jTjVj (k(6pA v +p8A v ) -SpS") ,(13) 
V-6u = 0, (14) 



V k 5u l 

SAij = [Si(kSi)j - AijAki) ,, ,, . (15) 



where 



We consider plane- wave disturbances (5u, dp) — exp (Xt + i£ ■ x)(a, a) prop- 
agating with wavevector £. In general a and a will be complex quantities. From 
the linearized equations, we obtain the following relations for A, a and a: 

Xai = Bia + Cija? , (16) 
i-a = 0, (17) 
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where 

Bi = J^2 (WAtj + i (kAij - Sij) e) , (18) 
dj = - (V jUi + »(£ ■ u)Sij) + (-fiM m £ k + iV k (jiMijki)) , (19) 

and 

M ijk i = (6i(j6i) k - AuAjk) . 

One can solve ( |l6| ) and (|l7|) for A for every wavevector £. The real part of A 
determines the growth or decay of the plane- wave disturbance with wavevector 
£, and the imaginary part determines the propagation of the disturbance. If the 
real part of A is positive for any £, we refer to this mode as linearly unstable, as 
this mode will grow rapidly in time. If, for a given background solution, ( p6| ) and 
( |l7| ) have unstable modes, then this solution will be termed linearly unstable. If, 
for a given background solution, there are unstable modes with arbitrarily short 
wavelengths, and if this background solution is a unique solution of the IPRF 
equations, then we will call these equations and the solution linearly ill-posed. 
Using the condition £ • a = we can eliminate a from the equation for A: 

Aa l = Dija j , (20) 
Da = (d ik -^jC kj , (21) 

a = -- S^V, (22) 
B ■ £ 

The eigenvalues of the matrix D (|2l]) determine the growth and propagation of 
the infinitesimal plane-wave disturbance. 

From ( pl| ) it can be shown that D has at least one zero eigenvalue Ai = 0. 
Since we are working in two dimensions, the remaining eigenvalue is equal to 
the trace of D: A2 = Tr(D). This trace is given by 

Tr(D) = Tr(C) - V _ , (23) 

= -i(u • + ^ (Vjt* + A lfe ^V% ; + Avb (A lk q k + ifiV k A lk )) . 
? • B 

where we have used the fact that the matrix is symmetric Ay = Auj) , is trace- 
free An = and further that Tr(Ai k A k j) = 1. 



We now consider the short wavelength (|£| — > 00) limit of (23). If we write 
X = then we can examine the powers of |£| on the right-hand side of (p3|). 

The leading order term in £ on the right-hand side goes as 0(|£| 2 ) and is real 
with coefficient 

(xiAi m Xm)(8ik ~ kA ik )x k Ajjx 3 
I 1 71 TTa \.. n .,fl 



{S pq - kA pq ) X Px 

This leading order term was considered by Schaeffer (p2|) in his analysis of (||) 
and (|^). The denominator is always positive for angles of friction ip > 0. The 
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numerator vanishes when \ points in the direction of the velocity characteristics 
(these lie at ±45° to the principal stress directions) or in the direction of the 
stress characteristics (these lie at angles of (±<p + 90°) /2 to the principal stress 
directions) of the background solution. To illustrate this, we have computed 
the eigenvalue A2 for a range of wavevectors for the approximate radial solution 
(given in the appendix) . The results are displayed in figure ^. In this contour 
plot the four wedges of unstable directions approach those given by ( p4| ) in the 
short wavelength limit. 

Plane wave disturbances with wavevectors £ that lie in directions between 
the stress and velocity characteristic directions will be unstable in the short 
wavelength limit as the real part of A2 will be positive in this limit. This short 
wavelength instability leads us to conclude that the two-dimensional granular 
flow equations are linearly ill-posed (this was first noted by Schaeffer [[j"2|). 



4 Singularly perturbed rigid-plastic flow equa- 
tions 

Attempts to study (||) and ([)]) numerically is problematic. Solutions of these 
equations obtained numerically are found to be grid-dependent due to the wedge 
of short wavelength instabilities ( J24|) found in the previous section. Here we pro- 
pose modifying the equations (|| -|9) in an effort to 'regularize' the instabilities 
( [l6| ) in some manner. We will consider the following equations: 

— + u^ju 1 = g l - V.j ( P S V - kpA v ) , (25) 
ViU* = eAp. (26) 

These equations are introduced as a tool to study shear-banding rather than 
an attempt to accurately model the physics of dense granular flows. However, 
it is worth noting some of the physical properties of these equations. These 
equations do not conserve mass as the density is held constant despite (^(|). 
They do satisfy a Drucker-Prager yield condition (Q) if one assumes a flow rule 
that is a natural generalization of (j^): 

(Tij - pSy = n (Vij - Tr{V)5 l3 /2) . (27) 

where fi is given by (||). Now, in this case 

fcV > {(Tij - pSrj) 2 - fe 2 p 2 (1 - <5 2 / 2), (28) 



where 5 = eAp/\\V\\. Note that equation (26) guarantees that \8\ < a/2- 

If 6 = 0, then the material is at yield and is incompressible, however as Ap 
increases the material relaxes from yield. If e is sufficiently small, one might 
expect that in some instances the solutions of (^5|) and ( |26| ) will approximate 
smooth solutions of (p6fjl7|). Nonetheless, we will see the inclusion of this term 



significantly alters the behavior of the linearized perturbation equations. 
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Figure 2: Contour plot of the real part of the eigenvalue A2 for the radial solution 
as a function of wavenumber (radial coordinate) and direction of the wave front 
(angular coordinate). Inside the wedges A2 is positive; outside the wedges, A2 
is negative. 
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The equations for plane-wave disturbances analogous to ( ^6[]l7j ) are: 

Xa t = li,<> ■ C,,<r. (29) 
£. a = -ie£ 2 a. (30) 



where B and C are given by ( |18| ) and ( |19| ) respectively. Using Q2S|) we can 
eliminate a directly from (|2^) to give: 

Xai = AjCtj = (Cij - i-^f)a 3 • (31) 

The leading order term on the right-hand side of ([5l]) in the short wavelength 
limit again goes as 0(|£| 2 ). This term is found to be 

M a = -j^2 M wx k x l - (32) 

The first step in analyzing the eigenvalues of D is to look at the eigenvalues 
of this leading order part M. In what follows we will set Fr — 1, and ignore 
the inertial terms which do not play a role in shear banding. If we write 

• eAp (33) 



then the trace of M is given by 

Tr(M) = -A(2 + 6^2 - S 2 cos(20)) . (34) 



and the determinant of M is 



A 2 



Det(M) = — (y 2 - S 2 cos(20) + 6) (35) 

where <f> is the angle of the wavevector £ to the principal stresses (here the stress 
is defined by equation (p7[)), and 

-Of 

Noting that \5\ < y/2, we see the trace of M is always negative. Further, 
Det(M) is non-negative (if \S\ > v2 — S 2 then Det(M) is strictly positive). 
Hence the matrix M has no positive eigenvalues, although one of the eigenvalues 
will vanish in the directions cos(20) = —S/y/2 — S 2 provided |5| < y2 — 8 2 . For 
a smooth background solution with |<5| <C 1, M will have a zero eigenvalue in 
four directions. 

In the short wavelength limit, the eigenvalues of D will tend to those of M 
(as these grow as |£| 2 ) except in the four directions where the determinant of 
M vanishes. Thus in this limit D will possess eigenvalues which are negative 
except very near these four isolated directions. This is in contrast to the earlier 
situation where equations Jl3| ) and (Q) possessed a wedge of unstable directions 
in the short wavelength limit. The possibility still remains that the equations 
and (pq) will be ill-posed near the four directions cos(2<^) = —5/y2 — S 2 . 
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5 Numerical solutions 



In this section we present numerical solutions of the equations ( J25| ) and (f26|). 
These equations are solved in a two-dimensional wedge-shaped hopper with fric- 
tional boundary conditions on the walls (6 = 6 W ) and stress-free inflow/outflow 
conditions on the arcs r = ro and r = r±. The radial solution described in 
the appendix is supplied as an initial condition. These computations have been 
using a finite element toolkit called Unstructured Grids (UG) which is freely 
distributed by the IWR at the University of Heidelberg (l) . 

The equations are discretized using a finite-volume scheme. Finite-volume 
methods have been successfully used to solve the Navier-Stokes equations in a 
variety of situations (Patankar pl[ ). We have chosen to use a control volume 
finite-element method due to Schneider and Raw ( |l5| ; see also Kariman and 
Schneider This method utilizes colocated primitive flow variables (i.e. aver- 
age stress and velocities which are discretized on the same grid) so is well suited 
to the formulation of the incompressible rigid-plastic flow equations presented 
in the previous section. The equations are discretized in time using a first-order 
backwards Euler scheme. 

The initial difficulty in solving these equations is the non-linearity in the 
stress terms on the right-hand side of ( psf ) . The convective terms in (|2^) are also 
non-linear. We will use a Newton linearization of these non-linear terms and 
solve the resulting equations using a Newton iteration. Within each Newton 
iteration, the resulting linear problem is solved using a biconjugate gradient 
method with multigrid preconditioning. 

The second difficulty is the growth of the instabilities and the development 
of shear bands. Once the instabilities develop a small time-step must be used. 
The shear-banding results in discontinuities developing across the wavefronts 
that grow from the instability. This is dealt with in an ad hoc manner; the 
grid is designed so that the discontinuities form parallel to element boundaries. 
The hopper has been discretized using grids of up to 4000 linear quadrangular 
elements. These fine grids are required to accurately resolve the right-hand side 
of (|2q) when e«l. On sufficiently fine grids, one finds that the wavelength of 
the instability is controlled by the magnitude of e. 

Figures || and [| depict the evolution of the flow from the smooth radial flow 
solution (see appendix). The figures show flows with (p = 34° and tp w = 9°. 
The parameter e is chosen to be 2 x 10~ 5 . The flows develop as follows: 

• Initially, plane wave-like wavefronts in the stress develop near the center of 
the hopper. This region corresponds to the location in the hopper where 
li is largest. 

• At these wavefronts, the tangential velocity component begins to develop 
jumps across the trough of each wavefront. 

• These instabilities develop further and further down the hopper, and the 
formation of jumps in the tangential velocity across the wavefronts follows. 
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• The growth of the wavefronts in amplitude does not increase without 
bound. Eventually, the growth appears to cut off. The wavefronts are 
observed to propagate down the hopper (see figure || where if = 34° and 
ip w = 9° but e is 1 x 10~ 4 so that the wavelength observed here is longer 
than in preceding figures). No steady state emerges. 

In the next section, we will analyze the linear stability of this system in further 
detail to try to understand these features. These features very much resemble 
shear bands, particularly the jumps in tangential velocity that occur across the 
wavefronts. It should be noted that the thickness across which this jump occurs 
seems to be grid-dependent. Thus, the shear bands are not fully resolved by the 
numerical method. 



6 Instabilities in the singularly perturbed sys- 
tem 

Including the singular perturbation in (pq), introduces a new scale, e, to the 
problem. The grouping ^/Jle gives a length-scale that appears in the modified 
eigenvalue problem (equation |3l| ). Indeed, the wavenumber of the instabilities 
in the hopper interior that are observed in the numerical studies above roughly 
satisfy /i|£| 2 ~ - 3> |£|, so this length-scale appears to be relevant for these 
features at least. This suggests that the eigenvalues of the following truncation 
of the full tangent matrix D should be sufficient to study the stability of these 
unstable modes: 

D' rj = My + - (kA ik - 6 ik ) X k Xr (36) 
e 

This truncation D' neglects the imaginary parts of D, and the contribution from 
the inertial terms. 

The trace and determinant of D' are given by: 

Tr(D') = Tr(M) - 1 ^1 + ^(y/2-6 2 cos(20) - 5)J , (37) 

Det(D') = Det(M) + -((2- S 2 ) cos 2 (20) + k^2 - 6 2 cos(20) + 5(6 - kj) . 

(38) 

In the short wavelength limit Tr(D') ~ Tr(M) and, unless Det(M) vanishes, 
Det(D') - Det(M). If, however, Det(M) vanishes then 

2A 

Det(D') 5(5 - k). (39) 

so that the determinant of D' will be negative if < 6 < k. In this case 
D' will possess a positive eigenvalue where Det(M) vanishes. Thus it seems 
plausible that under certain conditions ( p5| ) and ( p6| ) can become ill-posed along 
the directions where Det(M) may vanish. We should be cautious, however, in 
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t = 0.04 s 



t= 0.05 s 



Figure 3: Contour plots of the average stress in a wedge-shaped hopper showing 
the development of instabilities. Snapshots have been taken at 0.01 s intervals. 
These instabilities lead to shear-banding. 



12 



t = 0.00s 



t = 0.02s 



-Q.0& 0-S - 



=> a. 

a 

-D.15 



0.2 -> 

-I / 



t = 0.04s 



t = 0.06s 



Figure 4: Radial slices (constant 9) of the average stress, and the tangential 
velocity component across the band. The normal velocity components (not 
shown) are continuous across the bands. Snapshots are taken at 0.02 s intervals. 
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Figure 5: Radial slices (constant 9) of the average stress. This shows the wave- 
fronts propagating down the hopper. Snapshots are taken at 0.025 s intervals. 
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drawing conclusions about the short wavelength behavior of D, from the short 
wavelength behavior of D'. 

Away from the short wavelength limit, examination of Det(D') reveals the 
possibility of other instabilities. The existence of these instabilities depends on 
the size and sign of S. There are three possibilities which are classified below: 

1. > S > — 2( - v ^ +1 ^ ) '■ Det(D') is non- negative for all wavenumbers and all 
directions. 

2. 5 > k or — 2 ^y| +1 ^ > <5 ; Det(D') can be negative for wavenumbers that 
satisfy 

< Ae < 2(fc 2 - 4<S(<5 - k))/S{6 - k). 

Thus Det(D') is non-negative in the short wavelength limit (the problem 
is linearly well posed). 

3. k > 5 < 0: Det(D') can be negative at any wavelength (i.e. this suggests 
the problem may be linearly ill-posed). 

The situation is illustrated schematically in figure |^. 

This provides a qualitative explanation of the development of shear bands 
observed in the numerical simulations detailed in the previous section. The 
initial smooth background solution will have small, negative values of 6, so 
the problem is initially well-posed everywhere but unstable to perturbations of 
wavenumber |£| ~ l/^JJIe. As an example, the determinant of the full tangent 
matrix D has been computed for the radial background solution near the center 
of the hopper (see appendix), with tp — 34°, cp w = 9° and a value of e = 2 x 10 -5 . 
Illustrated in figure 0, we can see that the unstable region disappears as |£| — > oo, 
which is consistent with the behavior of D'. 

After a short time, an instability does develop with wavenumber |£| ~ 1/ ^/jje 
and wavevector at approximately 45° to the directions of principal stress. Once 
the disturbance grows sufficiently large, equations ( p9| - ^0| ) will cease to describe 
the behavior of the disturbance. Nonetheless, we can continue a qualitative 
description of the process using ( p9||3C| ) as a guide. As this disturbance grows 
then, S becomes a rapidly changing function of position in the hopper. In the 
troughs of the plane- wave- like pressure disturbance, S is positive, and the linear 
stability analysis would suggest that the solution becomes ill-posed in this case. 
Comparing this with the numerical simulations suggests that this ill-posedness 
in the troughs of the pressure disturbance causes a very strong gradient in the 
velocities to develop, which in turn drives S to zero. This feedback may provide 
a mechanism to regularize the ill-posedness. This needs to be studied further, 
by looking at length scales other than ^/Jle and by utilizing a non-linear analysis 
more suited to dealing with the inhomogeneities that develop as the instabilities 
grow. 
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Figure 6: Schematic stability phase diagram in terms of the Ae (which depends 
on frequency) and the parameter 5. There is now dependence on the direction of 
the wavevector is this plot; for points on the plot denoted unstable or ill-posed, 
there exists at least one direction with this property. 
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Figure 7: Contour plot of the real part of the eigenvalue A2 for the radial solution 
as a function of wavenumber (radial coordinate) and direction of the wave front 
(angular coordinate). The contour depicted is the A2 — contour. Inside the 
contour A2 is positive; outside this contour, A2 is negative. Here we can see the 
unstable region disappears as |£| — > 00. 
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7 Discussion 



We have presented new equations which exhibit shear-banding behavior. These 
equations are related to the IRPF equations by a singular perturbation. This 
singular perturbation regularizes the linear ill-posedness inherent in the IRPF 
equations by collapsing the wedge of wave directions that are unstable in the 
short wavelength limit. These resulting equations still possess instabilities, re- 
lated to those of the IPRF as seen in the linear stability analysis, and in nu- 
merical simulation these instabilities lead to shear-banding. The linear stability 
analysis suggests that the inhomogeneity that develops with these instabilities 
leads to further short wavelength instabilities, which cause these shear-band like 
features to develop. 

The wavelength of the instabilities observed in the singularly perturbed equa- 
tions was controlled by the time-scale e. The thickness across the jumps in tan- 
gential velocity across the wavefronts was found to be grid-dependent however. 
The linear stability analysis considered here was concerned with length scales 
^JJlI and did not predict a thickness across these jumps. However the analysis 
suggests that they are associated with the ill-posedness that can develop under 
certain circumstances. It is possible that the thickness of the shear bands here 
is zero; in this case, it would be impossible to fully resolve the band numeri- 
cally. Indeed, we were not able to resolve the bands in our numerical simulation. 
However, the linear analysis does suggest a mechanism whereby the ill-posedness 
may be shut-off as large velocity gradients appear in ill-posed regions. 

The stability properties of these equations is similar to that of the critical 
state plasticity studied by Schaeffer M| . Indeed, the parameter 8 (equation [53j) 
here plays a role similar to the dilation angle in CST. However, the equations 
proposed here are considerably simpler than the CST equations and may be a 
simpler system in which to study shear banding (the CST equations, however, 
will certainly provide a better model for granular flows). For instance, it may 
be possible to study the non-linear stability of this system to obtain a deeper 
understanding of the development of shear bands and their relationship to ill- 
posed behaviour. 
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A Approximate radial solutions 

For wedge-shaped hoppers Brennen and Pearce |2| and Kaza and Jackson Q, 
have examined steady-state plane-strain flows by expanding the velocities and 
the pressure in powers of (8/8 w ) where 6 W is the wall-angle of the wedge-shaped 
hopper. The steady-state equations are then solved order by order using a 
radial- flow ansatz. The first order terms in these expansions are similar to the 
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radial solutions of Jcnike for the quasistatic (Fr <C 1) steady-state equations 
(||) and (|9|). These approximations contain corrections due to the inclusion of 
the inertial terms when Fr ~ 1 . Working in polar coordinates (r, 9 where the 
hopper walls lie at ±9 W ) in the plane and using the Sokolovski parameterization 
of the stresses (p6|), the ansatz goes as follows: 



(40) 
(41) 



P = Po{r) + 0[—2 

\t>w 

* - *'(£) +0 (^<' 

where *f? is related to the stresses by 

* = 



2a 



r() 



a rr — aee 

The leading order terms in 9 are then given by 

F 2 1 



Po(r) = 



1 — sin tp 



w — 1 



r/r 
w — 1 



ro 



~ ) (f) 

(43) 



where *f? w is the value of \& at the wall (determined by the boundary condition 
(0)), 

- 2sinv? (i + «WM- 



(1 — sin ip) 



and 



^ 2 = 







\w — 


1) 



1- (10 

n 



w-l 1 



ui+2 



Here r = r$ and r = n are the arcs along which the average stress vanishes (the 
outflow and inflow respectively) . The corresponding terms in the velocity fields 
are also radial (to order (9/9 w ) ): 



fgr^F (\ 



ro 
r 



(44) 

v = 0, (45) 
where u and v are the radial and angular components of the velocity respectively. 
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